;-------------------------------------------------------------------------
; NCL script to read radon flux data
; Kai Zhang 2008-2010 
; http://en.wikipedia.org/wiki/ESRI_grid
;-------------------------------------------------------------------------

   load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"  

;;-------------------------------------------------------------------------
;; description 
;; 
;; E0: WCRP1995 emission 
;; E1: WCRP1995 magnitude with Schery1998 time information  
;; E2: Schery1998  Schery et al. (1998) scaled by 1.6 
;; E3: merged emission (china,europe,russia,usa,aus,schery1998) 
;; E10: Zhuo et al. (2008) + WCRP1995
;; E11: Zhuo et al. (2008) + WCRP1995 + Schery(2004)  
;; E12: Zhuo et al. (2008) + Schery1998 + Schery(2004)   
;; E13: 


;;-------------------------------------------------------------------------
;; default settings 
;;-------------------------------------------------------------------------
   IEMIS = 0
   land_limit = 0.500
   resolution = "63"
   land_limit_min = 0.1
   ntime = 12 

   r_min_flux = 0.001 ;; atoms/m2/s 
   r_max_flux = 2.000 ;; atoms/m2/s 
   l_min_flux = False ;; default 
   l_max_flux = False ;; default 

   l_landfrc_scale = False ;; default 
   l_WCRP_1995 = True  ;; default 

   l_new_russia = True
   l_new_usa    = True
   l_new_aus    = False 
   l_new_china  = True
   l_new_europe = True
   l_new_ocean_schery  = False ;; default 
   l_dflt_schery98 = False ;; default  

   l_sw1998ocn  = False 
   l_minfluxocn = False 

   l_zeropolarland = False 

   datapath = "./input/" 
   datapato = "./output/" 

;;-------------------------------------------------------------------------
;; user settings 
;;-------------------------------------------------------------------------
   load "./fr_option.ncl"

;;-------------------------------------------------------------------------
;; settings for each emission scheme 
;;-------------------------------------------------------------------------

   if(IEMIS.eq.0) then  
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s 
      r_max_flux = 5.000 ;; atoms/m2/s 
      land_limit = 0.500
      l_landfrc_scale = False ;; default 
      l_WCRP_1995 = True  ;; default  
      l_new_russia = False ;;True
      l_new_usa    = False ;;True
      l_new_aus    = False
      l_new_china  = False ;;True
      l_new_europe = False ;;True
      l_new_ocean_schery  = False ;; default 
      l_minfluxocn = True 
   end if 

   if(IEMIS.eq.1) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False ;; default
      l_WCRP_1995 = False ;; default
      l_dflt_schery98 = False 
      l_schery98_time = True
      l_new_russia = False ;;True
      l_new_usa    = False ;;True
      l_new_aus    = False
      l_new_china  = False ;;True
      l_new_europe = False ;;True
      l_new_ocean_schery  = False ;; default
      l_minfluxocn = True 
   end if

   if(IEMIS.eq.2) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False ;; default
      l_WCRP_1995 = False ;; default
      l_dflt_schery98 = True  
      l_new_russia = False ;;True
      l_new_usa    = False ;;True
      l_new_aus    = False
      l_new_china  = False ;;True
      l_new_europe = False ;;True
      l_new_ocean_schery  = False ;; default
      l_minfluxocn = True 
   end if

   if(IEMIS.eq.3) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False ;; default
      l_WCRP_1995 = False ;; default
      l_dflt_schery98 = True
      l_new_russia = True ;;False ;;True
      l_new_usa    = True ;;False ;;True
      l_new_aus    = True ;;False
      l_new_china  = True ;;False ;;True
      l_new_europe = True ;;False ;;True
      l_new_ocean_schery  = True ;;False ;; default
   end if

   if(IEMIS.eq.4) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False ;; default
      l_WCRP_1995 = False ;; default
      l_dflt_schery98 = True
      l_new_russia = False ;;True
      l_new_usa    = False ;;True
      l_new_aus    = False
      l_new_china  = False ;;True
      l_new_europe = False ;;True
      l_new_ocean_schery  = False ;; default
      l_sw1998ocn  = True 
   end if

   if(IEMIS.eq.5) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False ;; default
      l_WCRP_1995 = False ;; default
      l_dflt_schery98 = True
      l_new_russia = True ;;False ;;True
      l_new_usa    = True ;;False ;;True
      l_new_aus    = True ;;False
      l_new_china  = True ;;False ;;True
      l_new_europe = True ;;False ;;True
      l_new_ocean_schery  = True ;;False ;; default
      l_zeropolarland = True ;;False ;; default
   end if

   if(IEMIS.eq.11) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s 
      r_max_flux = 5.000 ;; atoms/m2/s 
      land_limit = 0.500
      l_landfrc_scale = False 
      l_new_ocean_schery  = True 
      l_WCRP_1995 = False  
   end if

   if(IEMIS.eq.21) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s 
      r_max_flux = 5.000 ;; atoms/m2/s 
      land_limit = 0.500
      l_landfrc_scale = False 
      l_new_ocean_schery  = True 
      l_WCRP_1995 = False  
      l_new_aus    = True
   end if

   if(IEMIS.eq.22) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s 
      r_max_flux = 5.000 ;; atoms/m2/s 
      land_limit = 0.500
      l_landfrc_scale = False 
      l_new_ocean_schery  = True 
      l_WCRP_1995 = False  
      l_new_aus    = True
      l_dflt_schery98 = True  
   end if

   if(IEMIS.eq.12) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.500
      l_landfrc_scale = False
      l_new_ocean_schery  = True
      l_WCRP_1995 = False  
      l_dflt_schery98 = True  
   end if

   if(IEMIS.eq.50) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.000 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.000 ;; not land limit  
      l_landfrc_scale = True 
      l_new_ocean_schery  = False ;; default
      l_WCRP_1995 = True  ;; default  
   end if

   if(IEMIS.eq.51) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.200
      l_landfrc_scale = True 
      l_new_ocean_schery  = True
      l_WCRP_1995 = False
   end if

   if(IEMIS.eq.52) then
      l_min_flux = True
      l_max_flux = True
      r_min_flux = 0.005 ;; atoms/m2/s
      r_max_flux = 5.000 ;; atoms/m2/s
      land_limit = 0.200
      l_landfrc_scale = True
      l_new_ocean_schery  = True
      l_WCRP_1995 = False
      l_dflt_schery98 = True
   end if


   if (l_min_flux) then
      print(" ")
      print(" smallest flux rate set to "+r_min_flux+" atoms/m2/s ")
      print(" ")
   end if

   print(" ")
   print(" set output resolution ")
   print(" ")

  ;pingo_grid = "r720x360.grid"
  ;gridtype = "05x05"

   pingo_grid = "t"+resolution+"grid"
   gridtype   = "T"+resolution

   print(" ")
   print(" gridtype: "+gridtype)
   print(" ")

   print(" ")
   print(" set file name ")
   print(" ")

   filanama = datapath+"emis."+"russia"+".nc"
   filanamb = datapath+"emis."+"europe"+".nc"
   filanamc = datapath+"emis."+"china" +".nc"
   filanamd = datapath+"emis."+"ocean" +".nc"
   filaname = datapath+"emis."+"schery"+".nc"
   filanamea= datapath+"emis."+"scherytime"+".nc"
   filanamf = datapath+"emis."+"usa"+".nc"
   filanamg = datapath+"emis."+"aus"+".nc"

   filenama = datapato+"out."+"russia"+"."+gridtype+".nc"
   filenamb = datapato+"out."+"europe"+"."+gridtype+".nc"
   filenamc = datapato+"out."+"china" +"."+gridtype+".nc"
   filenamd = datapato+"out."+"ocean" +"."+gridtype+".nc"
   filename = datapato+"out."+"schery"+"."+gridtype+".nc"
   filenamea= datapato+"out."+"scherytime"+"."+gridtype+".nc"
   filenamf = datapato+"out."+"usa"+"."+gridtype+"."+year+".nc"
   filenamg = datapato+"out."+"aus"+"."+gridtype+".nc"

   print(" output filenames: "+filenama)
   print("                   "+filenamb)
   print("                   "+filenamc)
   print("                   "+filenamd)
   print("                   "+filename)
   print("                   "+filenamea)
   print("                   "+filenamf)
   print("                   "+filenamg)

   print(" ") 
   print(" clean old files ")
   print(" ") 

   system(" rm " + filenama)
   system(" rm " + filenamb)
   system(" rm " + filenamc)
   system(" rm " + filenamd)
   system(" rm " + filename)
   system(" rm " + filenamea)
   system(" rm " + filenamf)
   system(" rm " + filenamg)

   print(" ")
   print(" remapping ")
   print(" ")

   print(" ")
   print(" ")
   system(" cdo remapcon,"+pingo_grid+"   "+filanama+ "  " + filenama) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamb+ "  " + filenamb) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamc+ "  " + filenamc) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamd+ "  " + filenamd) 
   system(" cdo remapcon,"+pingo_grid+"   "+filaname+ "  " + filename) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamea+"  " + filenamea) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamf+ "  " + filenamf) 
   system(" cdo remapcon,"+pingo_grid+"   "+filanamg+ "  " + filenamg) 
   print(" ")
   print(" ")

   filenamx = datapath+"T"+resolution+"_jan_surf.nc"

   print(" ")
   print(" read files ")
   print(" ")

   fila = addfile(filenama,"r") 
   filb = addfile(filenamb,"r") 
   filc = addfile(filenamc,"r") 
   fild = addfile(filenamd,"r") 
   filE = addfile(filename,"r") 
   filEa= addfile(filenamea,"r") 
   filf = addfile(filenamf,"r") 
   filg = addfile(filenamg,"r") 
   filx = addfile(filenamx,"r") 

   filenamo = datapato+"rnemis.E"+IEMIS+"."+gridtype+"."+year+".nc"

   system("rm "+filenamo)

   filo = addfile(filenamo,"c")

   print(" ")
   print(" read variables ")
   print(" ")

   vara = fila->rnemis 
   varb = filb->rnemis
   varc = filc->rnemis
   vard = fild->rnemis
   vare = filE->rnemis
   varea= filEa->rnemis
   varf = filf->rnemis
   varg = filg->rnemis
   oro  = filx->SLM

   lat = fila->lat
   lon = fila->lon

   nx = dimsizes(lat)
   ny = dimsizes(lon)

   print(" ")
   print(" define new vars ")
   print(" ")

;;-------------------------------------------------------------------------
;; define var 
;;-------------------------------------------------------------------------
   vary = new((/12,nx,ny/),"float")

   vary!0        = "time"
   vary&time      = fspan(1.,12.,12)
   vary&time@long_name = "time slice since Jan"

   vary!1        = "lat"
   vary&lat      = lat
   vary&lat@long_name = "latitude"

   vary!2        = "lon"
   vary&lon      = lon
   vary&lon@long_name = "longitude"

   vary@long_name = "radon flux"
   vary@units = "atmos/cm2/s"

   vary = -9999.

   print(" ")
   print(" combine emission over different regions ")
   print(" ")

   varz = vary(:,:,:)

   ;;printVarSummary(vara) 
   ;;printVarSummary(varz) 

   if (.not.l_new_ocean_schery) then 
      vard = (/ r_min_flux /)  
   end if 

;;-------------------------------------------------------------------------
;; set radon flux region by region 
;; ocean -> russia -> europe -> china -> others 
;;-------------------------------------------------------------------------

   if(l_new_ocean_schery) then
   do it = 0, ntime-1 
      varz(it,:,:) = where(.not.ismissing(vard(it,:,:)).and.oro(:,:).lt.0.5,vard(it,:,:),-9999.)  ;; ocean flux  
   end do 
   end if 

   if (l_new_russia) then 
      varz = where(.not.ismissing(vara),vara,varz)    ;; russia  
   end if 

   if (l_new_europe) then 
      varz = where(.not.ismissing(varb),varb,varz)    ;; europe 
   end if 

   if (l_new_usa) then 
      varz = where(.not.ismissing(varf),varf,varz)    ;; russia  
   end if 

   if (l_new_aus) then 
      varz = where(.not.ismissing(varg),varg,varz)    ;; aus
   end if 

   if (l_new_china) then 
      varz = where(.not.ismissing(varc),varc,varz)    ;; china 
   end if 


;;-------------------------------------------------------------------------
;; set other regions to Schery_1998 values 
;;-------------------------------------------------------------------------

   if (l_dflt_schery98) then 
      vare = vare / 21. / 1.6 
      print("                               ") 
      print(" scale Schery_1998 data by 1.6 ") 
      print(" ----------------------------- ") 
      print("                               ") 
      varz = where(varz.lt.0.,vare,varz)   
   end if 


;;-------------------------------------------------------------------------
;; reset radon flux if WCRP95 flux 
;;-------------------------------------------------------------------------

   if (l_WCRP_1995) then 
      print(" ")
      print(" pre-setup for WCRP95 emission")
      print(" ")
      varz = -9999.
   end if


;;-------------------------------------------------------------------------
;; scale factor for landfrac  
;;-------------------------------------------------------------------------

   scalef = dble2flt(oro) 

   print(" ")
   print(" set WCRP95 emission for other regions ")
   print(" ")


;;-------------------------------------------------------------------------
;; set radon flux for other regions  
;;-------------------------------------------------------------------------

   do im = 0, ntime-1
   do ix = 0, nx-1
   do iy = 0, ny-1

      ;; if emission flux not yet set 

      if (varz(im,ix,iy).lt.0.) then

         varz(im,ix,iy) =  (/ vard(im,ix,iy) /)  

         if (l_min_flux) then 
            varz(im,ix,iy) = r_min_flux 
         end if 

         ;; standard WCRP_1995 settings 
 
         if(lat(ix).ge.-60. .and. lat(ix).lt. 60. ) then 
            varz(im,ix,iy) = 1.
         end if

         if(lat(ix).ge.60. .and. lat(ix).lt. 70. ) then 
            varz(im,ix,iy) = 0.5
            if(lon(iy).ge.300. .and. lon(iy).lt. 360. ) then 
               varz(im,ix,iy) = 0.005 
            end if
         end if

         ;; set ocean flux 
         ;; zero flux when land fraction is too small 

         if(oro(ix,iy).lt.land_limit_min) then
            varz(im,ix,iy) = (/ vard(im,ix,iy) /)  
         end if


      end if

      ;; if l_landfrc_scale = T then scale flux according to land fraction, 

      if (l_landfrc_scale) then 
        varz(im,ix,iy) = varz(im,ix,iy) * scalef(ix,iy)  
      else   ;; l_landfrc_scale 
      end if ;; l_landfrc_scale  


      ;; reset emission flux over ocean when ORO less than land_limit  
      ;; this means not to use scaled land flux   
      ;; note that land_limit < land_limit_min

      if(oro(ix,iy).lt.land_limit) then
         if(l_minfluxocn) then 
            varz(im,ix,iy) = (/ r_min_flux /) 
         end if
         if(l_new_ocean_schery) then 
            varz(im,ix,iy) = (/ vard(im,ix,iy) /)  
         end if
         if(l_sw1998ocn) then 
            varz(im,ix,iy) = (/ vare(im,ix,iy) /)  
         end if
      end if

      if(oro(ix,iy).gt.land_limit) then
      if(l_zeropolarland) then  
         if(lat(ix).le.-60.) then
            if(im.ge.2.and.im.le.10) then 
               varz(im,ix,iy) = 0. 
            end if
         end if
      end if 
      end if 

   end do
   end do
   end do

;;-------------------------------------------------------------------------
;; recheck whether missing values and less-than-zero values exist, if so, zero flux at those locations  
;;-------------------------------------------------------------------------

   varz = where(ismissing(varz),0.,varz) 
   varz = where(varz.lt.0.,0.,varz) 


;;-------------------------------------------------------------------------
;; copy data to output variable 
;;-------------------------------------------------------------------------

   do im = 0, 11
      vary(im,:,:) = varz(im,:,:)
   end do

   if(IEMIS.eq.1) then
      vary = vary * varea 
   end if 

   print(" ")
   print(" output files ")
   print(" ")

   filo->rnemis=vary

   print(" ")
   print(" ")
   print(" ")
   print(" --- job successfully done --- ")
   print(" ")
   print(" ")
   print(" ")










